A  fast  solver  is  presented  for  the  solution  of  scattering  problems  in  which  the  scatterer  is 
a  relatively  thin,  elongated  object.  The  scheme  presented  here  is  a  version  of  an  algorithm 
previously  published  by  the  authors,  and  is  based  on  the  observation  that  under  certain 
conditions  and  with  certain  modifications,  the  scheme  will  retain  its  0(n )  CPU  time  estimate 
independently  of  the  size  of  the  scatterer  in  wavelengths.  The  performance  of  the  scheme 
is  illustrated  with  numerical  examples. 
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A  fast  direct  solver  for  scattering  problems 
involving  elongated  structures 

P.G.  Martinsson  and  V.  Rokhlin 


1.  Introduction 

A  standard  approach  to  the  numerical  solution  of  large-scale  scattering 
problems  is  to  reduce  the  problem  to  a  set  of  integral  equations  on  the 
boundary  of  the  scatterer  (or  scatterers),  discretize  the  integral  equations 
via  an  appropriate  numerical  scheme  (usually,  Nystrom  or  Galerkin)  and  to 
solve  the  resulting  system  of  linear  algebraic  equations  iteratively;  GMRES 
appears  to  be  the  preferred  iterative  solver  in  this  environment.  While  in 
many  situations  the  existing  methods  are  quite  satisfactory,  the  iterative 
methods  tend  to  display  slow  convergence  whenever  the  condition  number 
of  the  underlying  scattering  problem  is  high  ( i.e .  when  the  system  is  close  to 
a  resonance) .  Furthermore,  when  the  problem  to  be  solved  involves  a  single 
scattering  structure  irradiated  from  many  directions  (a  frequently  encoun¬ 
tered  situation),  iterative  methods  obtain  only  limited  advantage  from  the 
resulting  multiple  right-hand  sides. 

In  this  paper,  we  present  an  algorithm  for  the  direct  inversion  of  inte¬ 
gral  operators  associated  with  scattering  problems  involving  thin,  elongated 
scatterers.  The  computational  cost  of  the  scheme  is  O(N),  where  N  is  the 
number  of  nodes  in  the  discretization  of  the  surface  of  the  scatterer  (which 
in  this  case  tends  to  be  proportional  to  the  length  of  the  scatterer  in  wave¬ 
lengths).  The  algorithm  of  this  paper  should  be  viewed  as  a  development  of 
the  direct  inversion  scheme  found  in  [12].  A  class  of  problems  very  close  to 
the  one  addressed  here  is  treated  in  [14],  where  an  0(N  log2  N)  algorithm 
is  presented;  [14]  is  in  turn  an  extension  of  [8], 

This  paper  is  structured  as  follows:  Section  2  defines  a  model  problem 
and  lists  some  known  techniques  for  discretizing  and  solving  such  problems. 
Section  3  contains  the  analytical  apparatus  we  use  in  the  construction  of 
the  numerical  scheme  of  this  paper,  and  in  Section  4,  we  describe  the  “fast” 
algorithm.  Section  5  illustrates  the  performance  of  the  scheme  via  several 
numerical  examples  involving  scatterers  that  are  up  to  three  thousand  wave¬ 
lengths  in  size.  Finally,  Section  6  contains  generalizations  and  conclusions. 


2.  Preliminaries 

In  this  section,  we  describe  a  scattering  problem  that  will  serve  as  our 
model,  and  reduce  the  model  problem  to  an  integral  equation  of  the  sec¬ 
ond  kind  on  the  boundary  of  the  scatterer.  We  also  summarize  some  of 
the  standard  techniques  for  the  numerical  solution  of  the  obtained  integral 
equation. 
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2.1.  Integral  equations  resulting  from  scattering  problems.  For  defini¬ 
tiveness,  we  will  be  discussing  the  solution  of  the  exterior  Dirichlet  problem. 
Specifically,  given  a  complex  wavenumber  k  (which  is  assumed  to  have  a 
non-negative  imaginary  part)  and  a  subset  f 1  of  R2  bounded  by  a  suffi¬ 
ciently  smooth  Jordan  curve  T,  we  seek  the  function  ttout  :  R2  \  Q  — >  C 
satisfying  the  Helmholtz  equation 

(2.1)  — Aitol,t  —  k  uol,t  —  0, 


subject  to  the  radiation  condition 

(2.2)  f--Uont(t-x) 

for  any  x  € 

(2.3) 


lim 

t—*  OO 


oikt 


=  c(a) 


such  that  1 1  x*  1 1  =  1,  and  to  the  boundary  condition 

Uin  +  Uout  —  0; 


on  T,  with  ulu  :  r  — »  C  a  continuous  function. 

In  order  to  solve  the  boundary  value  problem  (2.1),  (2.2),  (2.3),  numeri¬ 
cally,  we  construct  the  Ansatz 

(2.4)  uout(x)  =  J  Ho(k\x  ~  2/1)  +  ikH0(k\x  -  y|)^j  a{y)  ds(y), 

where  n(y)  is  the  outward  pointing  unit  normal  of  F  at  the  point  y,  a  is  a 
function  to  be  determined,  and  Ho  is  the  Hankel  function  of  the  first  kind  of 
order  zero.  (This  choice  of  an  Ansatz,  which  leads  to  what  is  called  the  “com¬ 
bined  field  integral  equation”,  is  discussed  in  detail,  inter  aha,  in  [7].)  The 
function  uout  defined  by  (2.4)  automatically  satisfies  the  Helmholtz  equation 
(2.1)  with  the  radiation  condition  (2.2),  and  the  boundary  condition  (2.3) 
is  satisfied  if  and  only  if  a  satisfies  the  equation 


(2.5)  -  2ia(x)  +  ^~-^H0(k \x  -  y\)  +  ikH0(k\x  -  y|)^  a{y)  ds{y) 

for  all  x  e  T.  In  effect,  the  process  described  in  the  preceding  paragraph 
converts  the  boundary  value  problem  (2.1),  (2.2),  (2.3)  into  the  integral 
equation  (2.5).  The  purpose  of  this  paper  is  to  solve  the  equation  (2.5) 
rapidly  in  the  case  where  the  scatterer  is  a  relatively  thin,  elongated  object. 


2.2.  Nystrom  discretization  of  integral  equations.  Having  formulated 
a  scattering  problem  as  an  integral  equation  of  the  form 

(2.6)  v(x)  +  J  K(x,y)v(y)  ds(y)  =  f(x),  for  x  G  T, 

(as  discussed  in  Section  2.1),  the  next  step  is  to  approximate  the  continuum 
equation  (2.6)  by  a  system  of  N  linear  algebraic  equations.  There  are  several 
ways  of  doing  so;  throughout  this  paper,  we  have  utilized  a  version  of  the 
Nystrom  method.  This  technique  has  been  described  in  detail  elsewhere 
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(see,  e.g.  [2,  7]),  and  here  we  simply  observe  that  the  end  result  of  the 
process  is  a  set  of  N  linear  algebraic  equations  of  the  form 

(2.7)  (I  +  A)v^  =  f{N), 

where  A  is  a  matrix  that  approximates  the  integral  operator.  It  is  important 
to  note  that  when  the  Nystrom  method  is  used,  the  entries  of  the  matrix  A 
are  in  a  certain  sense  samples  of  the.  kernel  function  K(x,y).  To  be  precise, 
if  the  contour  is  discretized  into  the  points  {xi}^,  and  if  we  let  axj  denote 
the  ij- th  entry  of  A,  then  there  exist  functions  ip  and  ip  such  that 

(2.8)  aij  =  (p(x-j)K (.t,;,  x,j)ip(xj),  when  |i  —  j\  >  m. 

In  equation  (2.8),  the  matrix  elements  in  a  diagonal  band  of  width  m  are 
excluded  since  the  kernel  function  is  typically  singular  near  the  diagonal  and 
the  quadrature  weights  must  then  be  adjusted  to  retain  high  order  accuracy, 
see,  for  example,  [Kapur Rokhlin]  (typically,  m  <  10). 

2.3.  Iterative  solvers  for  the  linear  systems  resulting  from  scatter¬ 
ing  problems.  The  system  of  linear  algebraic  equations  (2.7)  is  commonly 
solved  using  an  iterative  solver  such  as  GMRES  in  conjunction  with  a  fast 
matrix- vector  multiplication  algorithm  such  as  the  Fast  Multipole  Method 
that  computes  a  matrix-vector  product  in  0(N)  operations.  If  the  iterative 
solver  needs  iVjter  iterations  to  converge  to  some  specified  accuracy,  then 
the  entire  solution  process  requires  0(iVjter  N)  floating  point  operations.  In 
many  environments,  N-lter  is  small  and  this  approach  works  very  well. 

The  principal  drawback  of  iterative  methods  in  the  context  of  high-frequency 
scattering  problems  is  that  such  problems  are  almost  always  close  to  being 
resonant,  in  which  case  a  large  number  of  iterations  may  be  needed  (since 
the  system  (2.7)  gets  ill-conditioned),  which  in  turn  slows  down  the  solution 
process.  Moreover,  it  is  common  that  a  sequence  of  scattering  problems  (cor¬ 
responding  to  a  collection  of  different  right-hand  sides)  needs  to  be  solved 
with  the  same  geometry,  which  is  a  situation  that  favors  direct  methods  over 
iterative  ones. 

2.4.  Direct  solvers  for  the  linear  systems  resulting  from  scattering 
problems.  If  Gaussian  elimination  or  a  similar  scheme  were  used  to  solve 
the  system  (2.7),  the  number  of  floating  point  operations  required  would  be 
0(iV'!),  which  would  severely  limit  the  size  of  problems  that  could  be  solved 
in  practice.  However,  several  schemes  for  the  acceleration  of  the  direct  so¬ 
lution  of  linear  systems  arising  from  integral  equations  have  been  proposed 
over  the  last  decade,  see  e.g.  [4],  [10],  [12].  Most  such  schemes  rely  on 
rank-deficiencies  in  the  off-diagonal  blocks  of  the  matrix,  which  means  that 
they  do  not  perform  asymptotically  fast  for  high-frequency  problems  (we 
say  that  a  method  is  “fast”  if  its  computational  speed  is  0(N  log9  N)  for 
some  integer  q).  An  algorithm  presented  in  [5]  reduces  the  computational 
complexity  to  0(N 3/2)  for  two-dimensional  Lippmann-Schwinger  problems 
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(and  would  yield  0(N2)  complexity  in  three  dimensions),  but  to  the  au¬ 
thors’  best  knowledge,  there  does  not  exist  a  “fast”  algorithm  (in  the  sense 
described  above)  for  the  direct  solution  of  high  frequency  problems  on  gen¬ 
eral  domains. 

However,  for  the  special  case  of  scattering  from  elongated  objects  (see 
Section  4  for  a  precise  definition),  it  was  demonstrated  in  [14]  that  the 
system  (2.7)  can  be  solved  directly  in  0(Nlog2  N)  operations  (an  earlier 
result  in  this  direction  was  given  in  [8]).  The  explanation  for  this  is  to  be 
found  in  the  observation  that  for  elongated  scatterers,  the  rank  of  interaction 
between  two  separated  pieces  of  the  scatterer  scales  as  the  logarithm  of  the 
size  of  the  pieces  (measured  in  wave  lengths),  rather  than  scaling  linearly 
with  the  size  of  the  pieces  as  would  be  the  case  for  a  general  scatterer. 
Versions  of  this  observation  are  formulated  in,  e.g.  [3],  [8],  [9],  and  [14], 

In  Section  3,  we  formalize  and  extend  some  of  the  observations  found  in 
[3],  [8],  [9],  [14],  The  result  will  enable  the  conversion  of  many  fast  inversion 
schemes  devised  for  non-oscillatory  problems  (such  as,  e.g.,  those  presented 
in  [4],  [10],  [12])  to  fast  algorithms  for  high-frequency  scattering  problems 
on  elongated  scatterers.  We  describe  such  a  modification  of  the  scheme  of 
[12]  in  Section  4,  and  illustrate  the  performance  of  the  resulting  scheme  on 
some  practical  examples  in  Section  5. 

2.5.  An  interpolation  result.  The  following  lemma  (that  will  be  used  in 
Section  3)  states  that  for  any  n-dimensional  linear  space  X  of  continuous 
functions  on  a  given  compact  set,  there  exists  a  set  of  n  interpolation  points 
such  that  any  function  in  X  can  be  interpolated  exactly  from  its  values  at 
those  n  points.  Moreover,  none  of  the  interpolation  coefficients  are  large. 
Lemma  1  is  a  special  case  of  Theorem  2  of  [13]. 

Lemma  1.  Suppose  that  Q  is  a  compact  set  in  M2,  J  is  a  positive  integer, 
and  that  xi are  continuous  complex-valued  funtions  on  Q.  Then 
there  exist  J  points  y\,. . .  ,yj  in  fl  and  J  functions  <p>\, ...  ,<pj  on  Q  such 
that,  for  j  —  1, ...  ,J  and  y  €  O, 

J 

(2-9)  Xjiv)  = 

<7=1 

and. 

(2.10)  \<Pj(y)\  <  i. 

3.  Low-rank  approximations 

In  this  section,  we  prove  certain  estimates  regarding  the  ranks  of  interac¬ 
tions  between  different  parts  of  an  elongated  scatterer.  The  starting  point 
is  Theorem  2  which  can  be  found  (in  a  somewhat  different  form)  in  [14].  It 
states  that  to  within  precision  e,  the  rank  of  interaction  between  the  two 
boxes  shown  in  Fig.  1,  of  lengths  L  and  fixed  widths,  is  0(log(fcL)  |  logs)2), 
as  L  — *  oo  and  e  — >  0.  A  proof  of  Theorem  2  can  be  found  in  Appendix  A 
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L  L 


Figure  1.  The  boxes  fig  and  fix  contain  source  and  target 
points,  respectively.  A  “typical”  source  point  y  and  a  “typ¬ 
ical”  target  point  x  are  marked  with  crosses.  The  point  y  is 
the  vertical  projection  of  y  onto  the  bottom  boundary.  The 
origin  in  the  coordinate  system  used  is  labeled  O. 

to  this  paper.  Theorems  3,  4,  5  are  the  purpose  of  this  section;  each  of  them 
can  be  viewed  as  a  modification  of  Theorem  2.  Theorems  3,  4,  5  will  be 
used  in  the  construction  of  low-rank  factorizations  of  off-diagonal  blocks  in 
the  matrices  to  be  inverted. 

In  the  following  theorem  (as  elsewhere  in  this  paper),  Ho  denotes  the 
Hankel  function  of  the  first  kind  of  order  zero. 

Theorem  2.  Suppose  that  L,  W ,  k,  e,  S  are  real  numbers  such  that  L  >  0, 
W  >  0,  k  >  0,  0  <  e  <  1/2,  S  >  0,  and  fis  =  [-L,  0]  x  [-W/2,  W/2], 
Ox  =  [5,  S  +  L]  x  [—W/2,  W/2\  are  two  boxes  in  R2  (see  Fig.  1). 

Then  there  exists  an  integer  J,  functions  ipj  :  fig  — >  C,  Xj  ■  fir  — >  C, 
with  j  =  1, 2,  ■  •  •  , and  a  real  number  c  that  does  not  depend  on  either  L 
or  £  such  that  if 

(3-1)  5>c|loge|/&, 

then 

./ 

(3.2)  H0(k\x  -  y\)  =  ^Ax)Xj{y)  +  E(x,  y),  for  all  x  <E  ftT,  y  €  fis, 

3= i 

where 

(3.3)  \E(x,y)  \  <  e. 

Furthermore,  there  exists  a  number  C  such  that,  as  L  — >  oo,  and  e  — >  0, 

(3.4)  J  <  Clog(A:L)| log | 2 . 

Finally,  the  functions  ifj  and  Xj  are  smooth  and  bounded  uniformly  in  L 
and  e. 

Remark  1.  In  Theorem  2,  there  is  a  requirement  that  the  boxes  must  be 
separated  by  a  distance  of  0{\  log  £\/k).  The  dependence  of  e  in  this  relation 
can  for  practical  purposes  be  ignored.  If  e  >  10~15,  then  a  separation 
S  —  1  / k  works  well.  The  reason  is  that  the  sum  (A.  13)  in  the  proof  converges 
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very  fast;  for  double  precision  accuracy,  no  more  than  15  terms  are  ever 
needed.  Another  consequence  of  this  fact  is  that  the  square  on  the  |loge| 
factor  in  (3.4)  is  generally  not  noticeable  in  practical  calculations. 

Remark  2.  Theorem  2  does  not  specify  how  the  interaction  rank  J  depends 
on  the  width  W  of  the  boxes.  In  fact,  it  grows  quite  rapidly  with  W,  and 
the  theorem  is  generally  not  practically  useful  for  W  larger  than  a  few  wave¬ 
lengths.  Some  actual  values  of  J  for  various  combinations  of  L,  W,  S,  and 
£,  are  given  in  Sec.  5.1. 

Theorem  3  sharpens  Theorem  2  by  stating  that  the  functions  tpj(x)  in  the 
expansion  (3.2)  can  in  fact  be  chosen  to  take  the  form  4>j(x)  =  Ho(k\x  - 
Dj\)  for  some  points  yj  €  Lls-  This  assertion  is  proved  by  applying  the 
interpolation  result  of  Lemma  1  to  the  expansion  (3.2). 

Theorem  3.  Let  L,  W,  k,  S,  e,  Lis,  and  J  be  as  in  Theorem  2.  Then 
there  exist  points  {y-j  }'-=l  in  Lis,  and  functions  {<Pj}j=i,  such  that, 

J 

(3.5)  H0(k\x~y\)  =  '52H0(k\x-yj\)<pj(y)  +  E(x,y)  for  x  £LlT,  y  e  Lls, 

3~  1 

where 

(3-6)  \E{x,y)\<e{J+l). 

Moreover,  for  j  —  1, . . . ,  J  and  y  e  Lls, 

(3.7)  \<Pj(y)\  <  1. 

Proof:  In  order  to  convert  the  results  of  Theorem  2  to  those  of  Theo¬ 
rem  3,  we  apply  the  interpolation  method  of  Lemma  1  to  the  set  of  func¬ 
tions  xi,  ■  ■  ■ ,  X.J-  This  yields  a  set  of  points  {yq}Jq=i  in  Lls,  and  interpolants 
{^,,}q=x  satisfying  (3.7).  Combining  (2.9)  and  (3.2),  we  find  that 

j  J 

(3.8)  H0{k\x  -  y\)  =  ^'(•x')Xj(y9)^9(%)  +  E(x,  y). 

1=1  '/=i 

Another  application  of  (3.2)  yields 

J 

(3.9)  H0{k\x  -  y|)  =  J2(H0(k\x  -  yq\)  -  E(x,yq))<pq(y)  +  E(x,y), 

q=l 

and  the  bound  (3.6)  is  seen  to  be  a  direct  consequence  of  (3.3)  and  (3.7).  □ 

The  following  theorem  is  the  principal  result  of  this  section.  It  is  simply  a 
corollary  of  Theorem  3  obtained  by  integrating  the  expansion  (3.5)  against 
a  charge  distribution. 

Theorem  4.  Let  L,  W ,  S,  k,  e,  Os ,  Llr r,  and  J  be  as  in  Theorem  2, 
and  let  the  functions  <p\, . . .  ,ipj  be  as  in  Theorem  3.  Furthermore,  for  any 
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function  a  £  Lx  (fis),  we  let  u  denote  the  potential  on  fix  induced  by  a  charge 
distribution  on  fis  of  density  a .  In  other  words,  for  x  £  fix, 

(3.10)  u(x)  =  /  H0(k\x-y\)cr(y)dA(y). 

■his 

Then,  if  we  define  for  j  —  1, . . . ,  J,  the  “equivalent”  charges 

(3.11)  Oj  =  f  <pj(y)cr(y)dA(y), 

■Jqs 

it  is  the  case  that,  for  x  £  fix, 

J 

(3.12)  u(x)  =  J^H0{k\x  -  yj\)aj  +  E(x), 

j  =  1 

where 

(3.13)  \E(x)\<(J+l)e  [  \a{y)\dA{y). 

Jns 

Proof:  Insert  the  expansion  (3.5)  into  (3.10).  □ 

Remark  3.  In  physical  terms,  Theorem  4  can  be  interpreted  as  follows: 
Suppose  that  fis  and  fix  are  as  in  Fig.  1.  Then  there  exist  J  points  inside  fis 
with  the  property  that  any  potential  in  fix  caused  by  a  charge  distribution  in 
fis ,  can  to  within  precision  e  be  replicated  by  placing  point  charges  at  those 
J  points  only.  Moreover,  the  number  of  points  depends  only  logarithmically 
on  the  lengths  of  the  boxes,  and  the  precision  required. 

The  functions  in  Theorem  4  form  a  basis  for  the  support  of  the 

integral  operator  that  maps  a  to  u  in  (3.10).  The  following  theorem  provides 
an  analogous  basis  for  the  range  of  the  same  integral  operator. 

Theorem  5.  Let  L,  W ,  S,  k,  e,  fis,  fix,  and  J  be  as  in  Theorem  f,  and  let 
for  any  a  £  L1  (fis),  the  potential  u  be  defined  by  (3.10).  Then  there  exist  J 
points  {xi}(_ j,  and  J  interpolation  functions  with  the  property  that, 

for  x  £  fix, 

J 

(3.14)  u{x)  —  ^2  Vi(x)u(xi)  +  E(x), 

i=  1 

where  E  satisfies  (3.13),  and,  for  i  =  1, . . . ,  J  and  x  £  fix, 

(3.15)  h(a)|  <  1. 

Proof:  We  let  an(l  {%}/=  i  be  defined  as  in  Theorem  4,  and  set 

J 

(3.16)  ue(x)  =  '22H0(k\x-yj\)crj. 

3= 1 

Thus  for  any  a,  the  function  u£  belongs  to  the  linear  span  of  the  functions 
{Ho(k\x  -  yj\)}j-i.  Lemma  1  then  directly  yields  the  desired  interpolation 
points  and  interpolation  functions.  □ 
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Remark  4.  Theorems  4  and  5  provide  bases  for  the  support  and  the  range  of 
the  integral  operator  in  (3.10),  respectively.  By  combining  the  two  theorems, 
we  obtain  a  factorization  of  the  integral  operator  into  three  factors  where  the 
middle  one  is  simply  the  restriction  of  the  original  operator  to  the  points 
{'Uj}j-i  and  {xi}Ji=l.  To  be  precise,  by  combining  (3.12)  and  (3.14),  we 
obtain  the  approximation 

(3.17)  u(x)  -  Ui(x)Ho(k\xi  -  yj\)  [  <pj{y)a{y)  dA(y)  +  E(x), 
ij= 1  JQs 

where  E  satisfies  (3.13). 

Remark  5.  While  Theorems  4  and  5  assert  the  existence  of  points  {yj}j= i 
and  {xi }/=1  such  that  (3.12)  and  (3.14)  hold,  they  do  not  say  anything  about 
how  to  actually  find  such  points.  Methods  for  this  task  are  described  in  [6] 
and  [13].  The  points  determined  by  these  methods  tend  to  lie  close  to  the 
boundaries  of  Os  and  Ox,  as  seen  in  Fig.  3(b).  In  fact,  it  is  possible  to  force 
the  points  to  lie  exactly  on  the  boundary,  although  one  must  then  use  both 
monopole  and  dipole  charges  (as  opposed  to  the  representation  (3.12)  that 
uses  monopoles  only). 

4.  Fast  direct  solvers  for  the  integral  equations  associated 
with  elongated  scatterers 

In  this  section,  we  describe  how  the  fast  solver  for  integral  equations 
presented  in  [12]  can  be  modified  in  order  to  obtain  a  “fast”  algorithm  for 
solving  the  integral  equation 

(4.1)  u(x)  +  J  H0{k\x  -  y\)u(y)ds(y)  =  f(x), 

where  /  is  a  given  function,  T  is  a  given  elongated  contour,  and  k  is  the 
wave-number.  For  simplicity,  we  assume  that  k  is  real  and  positive,  even 
though  the  scheme  works  for  complex  k  as  well,  as  long  as  the  imaginary 
part  of  k  is  non-negative.  The  discussion  is  framed  around  the  equation 

(4.1)  for  simplicity,  but  all  techniques  presented  apply  more  generally  to  the 
integral  equations  discussed  in  Section  2.1.  We  note  that  when  discussing 
the  asymptotic  cost  of  a  numerical  algorithm  for  high-frequency  problems, 
it  is  not  reasonable  to  fix  a  geometry  and  a  wave  number  and  let  the  number 
of  degrees  of  freedom  N  tend  to  infinity  independently.  Instead,  we  must  let 
the  size  of  the  object  (measured  in  wave  lengths  tend  to  infinity  in  proportion 
to  N. 

If  the  algorithm  of  [12]  were  to  be  applied  to  the  inversion  of  the  integral 
operator  in  (4.1)  unmodified,  its  asymptotic  cost  S  would  satisfy 

(4.2)  S  ~  c\N  +  C2(kL)3  as  iV, /cL —>  oo, 

where  L  is  the  diameter  of  F,  and  c\  and  C2  are  numbers  that  do  not  depend 
on  N,  k,  or  L.  It  turns  out  that  C2  is  quite  small,  so  for  scatterers  that 
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Figure  2.  The  contour  r  (drawn  with  a  solid  line)  is  ad¬ 
missible  provided  that:  (1)  the  width  W  of  the  enclosing 
rectangle  (drawn  with  a  dashed  line)  satisfies  W  <l/k,  and, 
(2),  there  is  a  number  no  such  if  m  denotes  the  number  of 
discretization  points  inside  a  box  of  length  D,  as  the  gray 
box  in  the  illustration,  then  m<nokD. 


are  at  most  a  few  hundred  wavelengths  large,  the  first  term  dominates  the 
computational  cost.  However,  the  second  term  eventually  will  dominate,  and 
will  limit  the  size  of  problems  that  can  be  considered.  In  this  section,  we 
demonstrate  how  the  scheme  can  be  modified  to  obtain  a  scheme  that  truly  is 
asymptotically  fast  in  the  sense  that  its  computational  cost  is  O(N)  —  0(kL) 
for  a  class  of  contours  satisfying  the  following  two  conditions: 

(1)  All  contours  in  the  class  fit  inside  rectangles  whose  shortest  side  is 
bounded  by  some  given  number  W.  (See  Remark  2,  and  Section  5.1, 
for  a  discussion  of  practical  values  of  W.) 

(2)  When  a  contour  in  the  class  is  discretized,  there  are  at  most  a  fixed 
number  no  of  points  per  wavelength,  counted  along  the  long  side  of 
the  rectangle  defined  in  (1). 

The  conditions  are  illustrated  in  Fig.  2.  We  will  from  this  point  on  assume 
that  the  long  side  of  the  box  is  oriented  along  the  horizontal  axis. 

Given  an  admissible  contour  T,  we  let  A  denote  the  N  x  N  matrix  that 
results  upon  Nystrom  discretization  of  the  integral  operator  in  equation 
(4.1).  If  the  points  in  the  discretization  are  ordered  based  on  their  horizontal 
coordinate  only,  it  follows  from  Theorem  4  that  the  rank  of  any  off-diagonal 
block  of  A  of  size  mxnis  0(log(m  +  n)).  This  shows  that  it  is  in  principle 
possible  for  the  algorithm  of  [12]  to  have  a  running  time  of  O(N).  To 
actually  achieve  this  speed,  we  must  furnish  the  algorithm  with  a.  technique 
for  rapidly  computing  bases  for  the  row  and  column  spaces  of  such  off- 
diagonal  matrices.  We  describe  such  a  technique  in  the  proof  of  the  following 
observation. 

Observation  6.  Let  T  denote  an  admissible  contour  of  horizontal  length 
L,  split  into  two  pieces  Ti  and  I^,  as  shown  in  Fig.  3(a).  Furthermore,  let 
m  and  n  denote  the  number  of  points  in  the  discretizations  of  Tj  and 
respectively,  let  B  denote  themxn  matrix  discretizing  the  integral  operator 
that  maps  a  charge  distribution  on  onto  a  potential  on  Fi,  and  let  {xj}^ 
denote  the  points  in  the  discretization  ofY\.  Then  for  any  given  precision 
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T2  Ti  T2 

(a) 


S  S 


Figure  3.  (a)  A  contour  r  =  Ti  +  T2-  Ti  is  drawn  with 
a  bold  line  and  T2  with  a  thin  line,  (b)  A  contour  T  = 
Ti  +  r close  +  Tfar-  Ti  is  drawn  with  a  bold  line,  rciOSe  with 
a  thin  line,  and  Ffar  with  a  dashed  line.  The  interaction 
between  Fj  and  rfar  is  replaced  by  the  interaction  between 
Ti  and  the  proxy  points  y^L\  y^R\  drawn  with  crosses. 


e,  there  exists  a  set  of  points  {s/y  }y=i  such  that  the  vectors 

(4.3) 


H0(k\xi  -  yi |) 

1 

g 

1 

_  H0(k\xm  —  2/1 1)  _ 

>  *  *  *  > 

_  H0(k\xm  -  yp  1)  _ 

constitute  a  basis  for  the  column  space  of  B,  accurate  to  within  precision  e. 
Moreover,  there  exists  a  number  C,  depending  on  neither  e,  nor  L,  such  that 

(4.4)  p  <  Clog(kL)  |  loge|2. 

The  cost  of  finding  the  points  {yjj^-i  is  0(p). 


Proof:  We  let  denote  the  points  provided  by  Theorem  4  in  order 

to  approximate  a  potential  associated  with  a  box  of  length  L  to  precision 
e,  and  let  S  be  the  corresponding  prescribed  separation  between  the  target 
box  and  the  source  box.  Given  this  separation  S,  we  let  rciose  denote  the 
part* of  T2  that  is  within  distance  S  of  Ti  in  the  horizontal  direction,  and 
let  Tfar  denote  the  rest,  see  Fig.  3(b).  Furthermore,  we  let  Bclose  denote  the 
matrix  mapping  a  charge  distribution  on  rciose  onto  a  potential  on  Fi,  and 
define  B^  analogously.  In  other  words,  the  matrices  B  and  [5close  -Bfar]  are 
identical  matrices  except  for  the  ordering  of  the  columns,  and 

(4.5)  Col(B)  =  Col(Bclose)  +  Col(Bfar), 

where  Col(B)  denotes  the  column  space  of  B.  We  let  {yjclose^}P^°so  de¬ 
note  the  points  in  the  discretization  of  rciose.  Then  the  vectors  {[Ho{k\xi  — 
2/]dOSe)|)]i=i}ji°r  form  a  basis  for  Bc iose.  It  remains  to  find  a  basis  of  the  re¬ 
quired  form  for  Col(Bfar).  But  this  is  exactly  what  Theorem  4  does.  Simply 
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£ 

=  io- 

-5 

£  = 

10"10 

kL: 

101 

~w 

103 

"IF" 

~w 

~W~ 

101 

To2- 

103 

~T(F~ 

~W 

~l¥~ 

kH  = 

1: 

15 

20 

24 

26 

28 

30 

29 

39 

48 

57 

64 

69 

kH  = 

2: 

20 

25 

28 

30 

32 

34 

36 

48 

57 

66 

73 

79 

kH  = 

3: 

24 

28 

33 

36 

38 

40 

44 

55 

64 

72 

80 

87 

Table  1.  The  table  displays  the  number  of  interpolation 
points  J  actually  required  for  the  interpolation  (3.12)  for  a 
separation  S  such  that  kS  —  1. 


enclose  the  right  part  of  Tfar  in  a  box  of  length  L,  and  let  de¬ 

note  the  associated  proxy  source  locations  obtained  from  the  set  {y]°^}/=1 
by  shifting  all  points  horizontally.  Likewise,  we  enclose  the  left  part  in  a 
similar  box  and  obtain  the  source  points  by  reflecting  the  points 

ab°ut  a  vertical  axis  and  then  shifting  them.  The  union  of  the 

three  sets  {yjclose)  }^°“° ,  {y^}j=1,  and  {yjR^}/=1  now  generates  a  basis  of 
the  required  form.  The  number  of  points  is  bounded  as  specified  by  (4.4) 
since  J  <  Clog(A:L)|  log£|2,  and  since  pclose  <  nokS  <  cno|  loge|.  □ 

5.  Numerical  examples 

5.1.  Numerical  estimation  of  interaction  ranks  between  elongated 
boxes.  Theorem  2  states  that  to  precision  £,  the  rank  of  interaction  J 
between  two  elongated  boxes  of  widths  W  and  lengths  L,  separated  by  a 
distance  S,  as  shown  in  Fig.  1,  satisfies 

(5-1)  J  <  C  log {kL)  (log(l/£))2, 

for  some  number  C,  as  L  — ►  oo  and  £  — >  0  (the  number  k  is  the  wave- 
number  in  the  Helmholtz  equation  (2.1)).  In  this  section,  we  present  the 
results  of  some  numerical  experiments  that  indicate  that  in  many  cases  of 
practical  interest,  J  is  quite  small,  typically  between  20  and  40.  Moreover, 
these  examples  indicate  that  even  though  Theorem  2  requires  the  separation 
distance  S  to  be  several  times  larger  than  the  wavelength,  a  separation  of 
only  one  half,  or  one  wavelength,  is  sufficient  for  five  or  ten  digits  of  relative 
accuracy. 

Table  1  gives  upper  bounds  for  the  interaction  ranks  for  a  number  of 
different  values  of  L  and  H  when  the  separation  distance  is  one  wave  length, 
while  Table  2  gives  the  corresponding  values  when  S  equals  one  half  wave 
length.  The  ranks  reported  in  the  tables  were  computed  using  the  methods 
described  in  [6]  and  [13]. 

5.2.  Example:  Scattering  off  an  elongated  closed  contour.  We  tested 
the  performance  of  the  inversion  scheme  described  in  Section  4  by  applying  it 
to  a  set  of  closed  contours  having  the  undulating  shape  shown  in  Fig.  3(a). 
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£ 

=  io- 

-5 

£  = 

t— 1 

o 

1 

O 

kL: 

101 

10y 

HF" 

"T(F 

"W 

~W~ 

101 

~W 

~w 

104 

To3- 

“lCF 

kH  —  1: 

20 

24 

29 

31 

33 

34 

35 

46 

55 

63 

70 

76 

kH  =  2: 

27 

33 

37 

40 

42 

44 

45 

55 

65 

73 

82 

87 

kH  =  3: 

35 

40 

44 

48 

49 

50 

55 

65 

75 

83 

92 

98 

Table  2.  The  table  displays  the  number  of  interpolation 
points  J  actually  required  for  the  interpolation  (3.12)  for  a 
separation  S  such  that  kS  —  0.5. 


^tot 

riwave 

4omp 

^apply 

ETes 

M 

800 

8 

2.38e0 

6.60e-3 

2.1e-6 

7.5e0 

1600 

16 

5.73e0 

1.54e-2 

3.3e-6 

1.7e0 

3200 

32 

1.27el 

3.26e-2 

3.5e-6 

3.6el 

6400 

64 

2.62el 

7.05e-2 

5.2e-6 

7.3el 

12800 

128 

5.39el 

1.41e-l 

7.5e-6 

1.5e2 

25600 

256 

1.19e2 

2.92e-l 

9.9e-6 

3.0e2 

51200 

512 

2.38e2 

5.89e-l 

1.2e-5 

5.9e2 

102400 

1024 

4.79e2 

1.19e-l 

1.2e-5 

1.2e3 

Table  3.  Results  of  the  numerical  experiment  described  in 
Section  5.2.  ntot  is  the  total  number  of  points  in  the  dis¬ 
cretization,  nwave  is  the  length  of  the  scatterer  in  wavelengths 
(its  width  is  1  for  all  experiments),  tcomp  is  the  time  in  sec¬ 
onds  required  to  compute  an  approximate  inverse,  tappiy  is 
the  time  in  seconds  required  to  apply  the  compressed  inverse 
to  a  vector,  Eies  is  the  residual  error  in  the  computed  inverse, 
and  M  is  the  amount  of  memory  (in  megabytes)  required  to 
compute  the  inverse. 


The  times  required  to  invert  the  integral  operator  in  (2.5)  on  a  desktop 
PC  with  a  3GHz  Pentium  IV  processor  and  3Gb  of  RAM  are  reported  in 
Table  3.  In  all  experiments,  the  vertical  size  of  the  contour  remained  fixed 
at  one  wave-length.  The  horizontal  size  varied  between  8  and  1024  wave¬ 
lengths  and  the  contour  was  discretized  using  50  points  per  wavelength.  The 
“geometric”  undulation  of  the  contour  has  precisely  twice  the  wave-length 
of  the  radiated  field. 

The  integral  equation  was  discretized  using  a  Nystrom  method  with  the 
6th  order  accurate  quadrature  rule  descriped  in  [11].  With  50  points  per 
wave-length,  this  yielded  a  discretization  error  of  10-10. 

5.3.  Example:  Scattering  off  a  one-dimensional  wavy  surface.  In 
this  section,  we  consider  the  problem  of  determining  the  scattered  wave 
from  an  infinite  contour  T  of  the  form  shown  in  Fig.  4.  The  contour  T 
coincides  with  the  x-axis  everywhere,  except  on  the  finite  piece  r+,  drawn 
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X  X  _ 

Figure  4.  Geometry  of  the  scattering  problem  described  in 
Section  5.3.  The  contour  F  is  drawn  with  a  solid  line;  the 
subset  T+  is  highlighted  using  a  bold  line.  The  reflection  T_ 
of  r+  is  drawn  with  a  dashed  line.  The  point  x  is  a  source, 
and  x _  its  mirror  image. 

with  a  bold  line  in  Fig.  4.  To  be  precise,  we  are  given  an  incoming  field  um 
of  the  form 

Ujn(:r)  =  Ho(k\x-x\), 

where  k  is  the  wave-number,  and  x  is  a  given  point  in  the  domain  Q  above 
T,  and  we  seek  a  field  uout  such  that 

(5-2)  -Auout  -  k2uoxl t  =  0,  on  Q, 

and 

(5-3)  wout  +  uin  =  0,  on  T, 

while  also  satisfying  the  radiation  condition  (2.2).  Letting  G(x,y)  denote 
the  kernel  of  the  integral  operator  in  (2.5),  we  make  the  Ansatz 

(5.4)  uout(x)  = 

[  G(x,y)cr(y)ds(y)  -  f  G{x,y~)  a(y)  ds(y)  -  H0(k\x  -  £_|), 
>/r+  J  r_ 

where  T_,  y_,  and  are  the  reflections  of  T+,  y,  and  x  in  the  2:- axis,  see 
Fig.  4.  With  this  Ansats,  the  Helmholtz’  equation  (5.2)  is  automatically 
satisfied,  and  the  boundary  condition  (5.3)  is  satisfied  on  r\r+.  To  ensure 
that  (5.3)  is  satisfied  on  r+  as  well,  it  is  necessary  that  a  satisfies  the 
equation 

(5.5)  -2ia(x)+[  G(x,y)  a(y)  ds(y)  -  [  G{x,yJ)  a(y)  ds(y)  - 

Jr+  J r_ 

Ho(k\x  —  x_|)  -  Ho(k\x  —  x|),  forxer+. 

The  times  required  to  invert  the  integral  operator  in  (2.5)  on  a  desktop 
PC  with  a  3GHz  Pentium  IV  processor  and  3Gb  of  RAM  using  the  method 
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ntot 

nwave 

tcomp 

^apply 

Eres 

M 

800 

11 

1.07e0 

2.20e-3 

9.2e-7 

2.6e0 

1600 

23 

2.32e0 

4.70e-3 

1.8e-7 

5.5e0 

3200 

46 

4.86e0 

9.60e-3 

7.7e-7 

l.lel 

6400 

91 

1.12el 

1.94e-2 

6.2e-7 

2.3el 

12800 

183 

2.30el 

3.90e-2 

8.8e-7 

4.7el 

25600 

366 

4.58el 

7.92e-2 

4.7e-6 

9.4el 

51200 

731 

1.01e2 

1.59e-l 

4.9e-6 

1.9e2 

102400 

1463 

2.02e2 

3.19e-l 

6.3e-6 

3.8e2 

204800 

2926 

4.06e2 

6.36e-l 

1.7e-5 

7.6e2 

Table  4.  Computational  timings  for  the  numerical  experi¬ 
ment  described  in  Section  5.3.  The  notation  is  the  same  as 
in  Table  3. 


described  in  Section  4  are  reported  in  Table  4.  In  all  experiments  reported, 
the  contour  T+  extends  for  half  a  wave-length  above  the  x-axis,  and  for  nwave 
wave-lengths  in  the  horizontal  direction,  where  10  <  nwave  <  3000.  The 
“geometric”  undulation  of  the  wave  pattern  has  about  the  same  wavelength 
as  the  radiating  field.  The  integral  equation  was  discretized  using  the  same 
techniques  that  were  described  in  Section  5.2. 

A  benefit  to  using  direct  methods  in  the  study  of  scattering  problems, 
is  that  one  frequently  is  interested  in  computing  the  scattered  field  for  a 
large  set  of  incoming  wave  patterns.  In  the  direct  method  presented  here, 
such  calculations  are  very  inexpensive  once  the  approximate  inverse  has  been 
calculated.  As  an  example,  we  computed  the  inverse  for  the  integral  operator 
corresponding  to  a  wavy  surface  such  as  the  one  depicted  in  Fig.  5,  but  600 
wave  lengths  across.  We  then  computed  the  reflected  field  for  a  set  of  sources 
located  at  various  positions  on  the  semi-circle  in  Fig.  5  and  computed  the 
energy  of  the  reflected  field  back  at  the  source  point.  This  resulted  in  a  graph 
showing  the  amount  of  energy  reflected  back  to  the  source  as  a  function  of 
location  on  the  semi-circle,  given  in  Fig.  6.  In  principle,  the  computation  of 
the  reflected  energy  for  each  source  location  requires  the  solution  of  a  dense 
linear  system  involving  40000  unknowns  (in  complex  arithmetics).  Since 
we  had  access  to  the  inverse  (which  took  about  a  minute  to  compute),  each 
computation  took  less  than  one  tenth  of  a  second  on  a  desktop  PC. 

6.  Conclusions 

In  this  paper,  we  have  demonstrated  that  the  integral  operators  associated 
with  acoustic  scattering  from  elongated  scatterers  in  two  dimensions  can 
be  inverted  very  rapidly.  Numerical  examples  show  that  for  a  scatterer 
3000  wavelengths  long,  the  forward  scattering  problem  can  be  solved  in 
minutes  on  a  desktop  PC  (at  five  digit  accuracy).  After  the  first  solve,  a 
compressed  representation  of  the  inverse  of  the  integral  operator  is  available, 
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Figure  5.  Computation  of  the  reflection  pattern  of  a  wavy 
surface.  A  radiating  source  is  placed  at  the  point  marked 
x,  the  reflection  problem  is  solved,  and  the  intensity  of  the 
reflected  wave  at  the  point  x  is  computed.  A  plot  of  this 
intensity,  as  a  function  of  the  angle  between  the  dashed  line 
and  the  positive  x-axis,  is  given  in  Fig.  6.  (Note  that  in  the 
experiment  reported,  the  wavy  surface  was  600  wave  lengths 
long  rather  than  the  11  shown  here.) 
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Figure  6.  This  plot  is  described  in  the  caption  of  Fig.  5. 
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and  additional  forward  scattering  problems  can  be  solved  in  less  than  one 
second  each. 

The  techniques  presented  can  with  minor  modifications  be  applied  to 
electromagnetic  scattering  problems  in  both  two  and  three  dimensions,  as 
long  as  the  scatterer  can  be  contained  in  a  long,  straight  cylinder  whose 
diameter  as  at  most  a  few  wave-lengths. 
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Appendix  A.  Proof  of  Theorem  2: 

The  proof  of  Theorem  2  given  in  this  appendix  relies  on  the  fact  that 
any  Hankel  function  Hn  can  be  approximated  by  a  sum  of  exponentials. 
The  number  of  terms  required  to  obtain  an  accuracy  of  £  on  the  interval 
[n,  n  +  L\  scales  as  log  L  |  logs).  This  result  is  stated  in  detail  in  Lemma  7. 
From  Lemma  7,  the  special  case  of  Theorem  2  in  which  the  two  boxes  JT2s  and 
f>T  have  zero  height  follows  immediately;  we  formulate  it  as  Lemma  8.  The 
proof  that  the  results  of  Lemma  8  can  be  generalized  to  the  environment  of 
Theorem  2  is  given  at  the  end  of  this  appendix. 

Lemma  7  is  a  slight  generalization  of  a  result  in  [8], 

Lemma  7.  Let  L  be  any  positive  number  larger  than  2,  let  n  be  any  integer, 
and  let  e  be  a  given  precision  such  that  0  <  e  <  1/2.  Then  there  exists 
an  integer  J,  a  number  C,  complex  numbers  {otj,/3j}j-i  o.nd  non-negative 
real  numbers  {tj}j=i  with  the  following  properties:  The  Hankel  function  Hn 
allows  the  expansion 

J 

(A.l)  Hn(r)  =  £  (ajeir  +  (3j)  e~^  +  E(r), 

j= l 

where 

(A. 2)  |i?(r)|  <  £,  forn<r<n  +  L. 

Moreover,  the  number  C  does  not  depend  on  either  L,  e,  or  n,  and 
(A. 3)  J  <  CTogL|log£|, 

Proof:  Formula  (9.1.22)  of  [1],  says  that 

(A.4)  Hn(r )  =  \E'„(r)  +  x„(r), 

where 

1  r°° 

(A. 5)  Vkn(r)  =  --  /  e~rsmht  (sm{mr)e-nt  +  i{ent  +  e~nt  cos(nTr)))  dt, 
n  Jo 
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and 

(A. 6)  Xn(r)  =  -  r  eirsine  e~ind  d6. 

Jo 

We  convert  the  integral  representation  of  to  an  expansion  of  the  form 
(A.l)  by  simply  estimating  the  integral  by  a  quadrature  rule  supported  on 
some  points  To  prove  that  the  number  of  points  J  required  in 

order  to  obtain  an  accuracy  of  e  is  bounded  as  prescribed  by  (A. 3),  we  note 
that,  (i),  the  integrand  in  the  representation  (A.5)  decays  exponentially 
(actually  super-exponentially),  and,  (ii),  the  derivative  of  the  integrand  at 
the  origin  is  bounded  by  the  condition  that  r  satisfy  1  <  r  <  1  +  L,  see  [15] 
or  [16]. 

It  remains  to  derive  an  expansion  for  Xn-  To  this  end,  we  first  effect  the 
change  of  variables  t  =  sin#  in  (A. 6),  whence 

(A.7)  Xn(r)  =  [  e'Ttxn(t)dt, 

Jo 

where 

(A.8)  Xn(t)  =  ^==|  ((V/T^2  +  it)»  +  +  it)n )  . 

Since  the  integrand  in  (A.7)  is  holomorphic  in  the  set  (0, 1)  x  (0,  oo),  we  can 
change  the  path  of  integration  from  the  horizontal  line  {t:0<t<l}to 
the  vertical  lines  {it :  0  <  t  <  oo},  and  (1  +  it  :  0  <  t  <  oo}.  Taking  into 
account  the  pole  at  t  =  1,  we  then  find  that,  for  some  number  c, 

r  oo  roc 

(A.9)  Xn(r)  =  /  e~rtxn{it)dt-  /  e*r_riXn(l  +  it)  dt  +  c e" . 

Jo  Jo 

The  integral  representation  (A.9)  can  now  be  converted  to  an  expansion  of 
the  form  (A.l)  via  the  quadrature  procedure  already  described  for  □ 

Lemma  8.  Let  L  be  any  positive  number  larger  than  2,  let  n  be  any  integer, 
let  s  be  a  given  precision  such  that  0  <  e  <  1/2,  and  consider  the  two  sets 
fis  =  {( yi ,  0)  :  -L  <  y\  <  0}  and  flj  =  {(^l,  0)  :  n  <  x\  <  n+L}.  Then, 
there  exists  an  integer  J,  a  number  C,  and  functions  {ijjj,  Xj}j= i  with  the 
following  properties:  The  Hankel  function  Hn  allows  the  expansion 

j 

(A.10)  Hn{\x-y\)  =  '^2ipj(x)xj(y)  +  E(x,y),  for  x  £  y  e  tts, 
j= l 

where 

(A. 11)  J  <  Clog L\ log e|, 

and 

(A-12)  \E(x,y)\  <  e. 

Moreover,  the  number  C  does  not  depend  on  either  L,  e,  or  n. 
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Proof:  In  the  geometry  given,  Hn{\x  —  y|)  =  Hn(x j  —  y\),  where  n  < 
x\  —  y\  <n  +  2L.  Thus,  the  formula  (A.l),  which  expresses  Hn  in  terms  of 
exponentials  only,  applies.  The  factorization  (A.  10)  follows  immediately.  □ 

Proof  of  Theorem  2:  Since  k  is  simply  a  scaling  parameter  for  the  geom¬ 
etry  of  the  problem,  we  can  for  the  purposes  of  the  proof  set  k  =  1. 

We  first  construct  an  expansion  of  the  form  (3.2)  that  is  valid  in  the 
special  case  where  x  lies  on  the  bottom  boundary  of  fix-  Picking  an  arbitrary 
y  =  (y i,  2/2)  £  Os ,  we  set  y  =  (y\ ,  —W/2),  as  shown  in  Fig.  1,  and  apply  the 
Graf  addition  theorem  (formula  (9.1.79)  of  [1])  to  obtain 

OO 

(A. 13)  H0(\x-y\)=  £  inHn(\x  —  y\)Jn{\y  —  y\). 


Since  \y  -  y\  <  W,  this  series  converges  very  fast;  to  be  precise,  we  can  for 
any  e  find  an  integer  N  such  that  N  ~  |  loge|,  and 


(A- 14) 


N 

Ho(\x  -  y\)  -  J2  inHn{\x-y\)Jn{\y-y\) 

n=-N 


<  £ 


for  any  y  €  Os,  and  any  x  on  the  bottom  of  fix-  Provided  that  S  >  N, 
Lemma  8  then  applies  to  each  term  Hn{ \x  —  yj)  in  (A. 14),  resulting  in  the 
desired  expansion.  The  condition  (3.1)  follows  since  N  ~  |loge|.  (See 
Remark  1  for  some  notes  on  the  truncation  of  the  sum  (A.  13).) 

Next  we  note  that  similar  expansions  can  be  constructed  for  x  on  all 
boundaries  of  fix  (the  top  boundary  is  of  course  entirely  analogous  to  the 
bottom  one,  and  the  corresponding  proof  for  the  vertical  boundaries  is  triv¬ 
ial).  Similarly,  it  can  be  shown  that  there  exist  expansions  of  the  desired 
form  for  the  normal  derivative  of  u  on  all  boundaries  of  Ox- 

Finally  we  note  that  when  we  know  the  value  of  u  and  its  normal  derivative 
on  the  boundary  of  Ox,  the  value  at  any  point  x  G  Ox  can  be  constructed 
via  Green’s  formula: 

(A.15) 

"(I)  =  5 L-r  “ zl)) "W  ~ HMx ~ 2|)*fe“(2))  ds(z) ' 

valid  for  any  function  u  satisfying  —A u  —  u  —  0  on  Ox-  □ 
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